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We study the Anderson localization of atomic gases exposed to three-dimensional optical speckles 
by analyzing the statistics of the energy-level spacings. This method allows us to consider realistic 
models of the speckle patterns, taking into account the strongly anisotropic correlations which are 
realized in concrete experimental conhgurations. We first compute the mobility edge Ec of a speckle 
pattern created using a single laser beam. We find that Ec drifts when we vary the anisotropy of the 
speckle grains, going from higher values when the speckles are squeezed along the beam propagation 
axis to lower values when they are elongated. We also consider the case where two speckle patterns 
are superimposed, forming interference fringes, and we find that Ec is increased compared to the case 
of idealized isotropic disorder. We discuss the important implications of our findings for cold-atom 
experiments. 

PACS numbers: 03.75.-b, 67.85.-d,05.60.Gg 


Anderson localization is the complete suppression of 
wave diffusion due to destructive interferences induced 
by sufficiently strong disorder [I|. It was first discussed 
by Anderson in 1958 Q and has been observed (only 
much lat^ in various physical systems, includiM light 
waves sound waves Q, and microwaves Q, and 

also in experiments performed with ultracold gases, 
first implementing an effective Anderson model Q, and 
then obser ving the localization of matter waves in one 
dimension |lOl ini and in three dimensions [3 M- 
Recently, transverse Anderson localization has been 
realized in randomized optical fibers 0, paving the 
way to potential applications in biological and medical 
imaging [3|- 

The key quantity which characterizes Anderson local¬ 
ization in three-dimensional quantum systems is the 
mobility edge Ec, which is the energy threshold that 
separates the localized states (with energy E < Ec) from 
the delocalized ones (with energy E > Ec) [3|- Many 
accurate theoretical predictions for the value of Ec exist, 
but most of them re gard simplified toy models defined 
on a discrete lattice |l7l . 0 • These lattice models do 
not describe the spatial correlations, and their possible 
anisotropy, of the disorder present in the physical 
systems where Anderson localization has been observed. 
In fact, these features are expected to have a profound 
impact on the Anderson transition. For example, it is 
known that due to finite spatial correlations an effective 
mobility edge exists also in low-dimensional systems [3- 
0, while for uncorrelated disorder all states would 
be localized [0. According to recent results [l^l, in 
continuous-space systems localization does not occur if 
the disorder correlation length vanishes, even for strong 
disorder. It is also known that the structure of the 
spatial correlations changes drastically the localization 
length and the transport properties [0 0 . 

The experiments performed with ultracold atoms are 
emerging as the ideal experimental setup to study 
Anderson localization [27l - l0 . Different from other 
condensed-matter systems, atomic gases are not affected 
by absorption effects and permit us to suppress the 


interactions. Furthermore, by shining coherent light 
through diffusive surfaces, experimentalists are able to 
create three-dimensional disordered profiles (typically 
referred to as optical speckle patterns) with tunable 
intensity and to manipulate the structure of their spatial 
correlations. 

In this Rapid Communication, we investigate the 
Anderson localization of noninteracting atomic gases 
moving in three-dimensional optical speckles. We 
determine the single-particle energy spectrum using 
large-scale diagonalization algorithms. Then, by per¬ 
forming a statistical analysis of the spacings between 
consecutive energy levels, we locate the mobility edge. 
The study of the level-spacing statistics lies at the 
heart of random-matrix and quantum-chaos theories. 
It has permitted us to interpret the complex spectra 
of large nuclei, atoms, and molecular systems [sol . [3l|. 
More recently, it has been employed in the analysis 
of the Google Matrix [32|> 101 • Quantum-chaos theory 
provides a universal basis-independent criterion for 
the localization transition. One has to identify two 
kinds of level-spacing distributions, namely, the Wigner- 
Dyson distribution characteristic of ergodic chaotic 
systems, and the Poisson distribution characteristic of 
localized quantum systems. This method has allowed 
researchers to locate the localization transition in 
noninteracting three-dimensional lattice models (both 
isotropic and anisotropic) [33 - 10 . and, more recently, 
also in interacting one-dimensional spin systems |38l - 
0. In the present study this criterion is used to 
investigate the Anderson localization of matter waves, 
setting the basis for future investigations of many-body 
localization in interacting three-dimensional Fermi gases. 

First, we consider the experimental configuration 
with a single speckle pattern created by shining a laser 
through a diffusive plate. In this case the spatial correla¬ 
tions of the disorder are intrinsically strongly anisotropic, 
with cylindrical symmetry around the beam propagation 
axis. We find that, when the speckles are elongated 
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FIG. 1: (color online) (a) Intensity profile of a speckle pattern 
measured on a plane orthogonal to the beam propagation axis 
z. (b) Elongated speckle pattern with anisotropy a^/cr = 6, 
measured along a plane containing the beam axis (c) Pro¬ 
file resulting from two orthogonally crossed speckle patterns 
(see text) measured on a plane containing the second principal 
axis, (d) Representation of the two speckle-patterns configu¬ 
ration, indicating the propagation directions of the first beam 
{z) and the second beam (y). The red arrows indicate the F* 
principal axis {z') and the 2“'^ p. a. (y'). The x-axis enters 
the sheet plane. 


along the axis, which is the typical experimental sit¬ 
uation, the mobility edge is only moderately reduced 
compared to the idealized models of disorder with 
spherically symmetric correlations. This unexpected 
result indicates that the experimental setup with a single 
speckle pattern is quite suitable to investigate Anderson 
localization, despite the strong disorder anisotropy. We 
also consider the case where two orthogonal speckle pat¬ 
terns are coherently superimposed. This setup, which 
was originally implemented to avoid the large axial 
correlation length of the single-pattern configuration, 
generates an intricate correlation structure, with rapid 
oscillations of the external field due to interference 
fringes (see Fig. [T|) [d^. In this case we find that the 
mobility edge is higher than for isotropic disorder, and 
is similarly to the case of a single speckle pattern with 
axially squeezed speckle grains. This means that the 
two-pattern configuration provides experimentalists a 
handle to shift upwards the position of the mobility edge. 

The first step in the determination of Ec is to com¬ 
pute the spectrum of the single-particle Hamiltonian 
H = — -|- H(r), where H is the reduced Planck’s 


constant, m is the atom’s mass, and H(r) is the dis¬ 
ordered potential experienced by the atoms exposed 
to optical speckle patterns. We consider a large box 
with periodic boundary conditions, which has a cubic 
shape (of size L) and parallelepiped shape, for isotropic 
and anisotropic speckles, respectively. We tackle this 
challenging computational task by representing H in 
momentum space, truncating the Fourier expansion 
at a large wave vector, carefully analyzing that the 
basis-truncation error is smaller than the final statistical 
uncertainty. To compute the eigenvalues we employ 
advanced numerical libraries for high-performance 
computers with shared-memory architectures [d^. For 
more details on the Hamiltonian representation and 
on the numerical diagonalization procedure, see the 
Supplemental Material M\. 

If the speckle field is blue detuned with respect to the 
atomic transitions, it generates a repulsive potential 
with an exponential probability distribution of the local 
intensity, which reads Phd{V) = exp {—V/Vq) /Vq, if 
the intensity is H > 0 and P{V) = 0 otherwise. Thus, 
the potential has the lower bound H(r) = 0, while it 
is unbounded from above. The disorder strength is 
determined by the energy scale Vq , which is equal to the 
spatial average of the potential, Vq = (^(r)) and also to 
its standard deviation, so that = (H(r)^) — (I/(r))^. 
For sufficiently large systems the disorder is self¬ 
averaging, and the spatial average coincides with 
the average over disorder realizations. Another fun¬ 
damental property which characterizes the speckle 
pattern is the two-point spatial correlation function 
F(r) = (H(r' -|-r)I/(r')) — 1. After averaging over 

the position of the first point r', it depends on the 
relative (vector) distance r. 

In order to make a direct comparison with a previous 
theoretical study based on transfer-matrix theory [i^ . 
we first consider an idealized isotropic model of the 
speckle pattern with a spherically symmetric correlation 
function that reads F’®°(r) = [sin(r/(T)/(r/cr)]^ (see 
inset in Fig. El). The parameter a fixes the length scale 
of the spatial correlations and therefore the typical 
grain size [43 |. An efficient numerical algorithm to 
generate isotropic speckle patterns is described in details 
in Refs. [13, [3. We determine the energy spectrum 
of a large number of realizations of the speckle pat¬ 
tern [15 . In the high-energy regime, the energy levels 
En (listed in ascending order) fluctuate, avoiding each 
other, signaling the level repulsion typical of delocalized 
chaotic systems. The distribution of the level spacings 
bn = En+i — En should Correspond to the statistics 
of random-matrix theory (in particular, to the Gaus¬ 
sian orthogonal ensemble), namely, the Wigner-Dyson 
distribution. Instead, in the low-energy regime the 
energy levels easily approach each other like independent 
random variables. This is a consequence of the localized 
character of the corresponding wave functions. In this 
regime the level spacings follow a Poisson distribution. 
In order to identify the two statistical distributions 
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FIG. 2: (color online) Two-point spatial correlation functions 
of the disorder. The continuous curves represent the analyt¬ 
ical formulas and the symbols represent the correlation 
measured on the speckle patterns generated numerically. The 
inset shows the correlation function of a single speckle-pattern 
along the beam axis 2 for elongated speckle grains (p^ja = 3) 
and squeezed speckle grains (pzjo = 1/3), radial correlation, 
and isotropic correlation of idealized spherically-symmetric 
speckle patterns. The main panel shows the correlation of 
crossed speckle patterns along the first ( 2 ') and the second 
principal axis (y') and along the orthogonal axis x. 



FIG. 3: (color online). The main panel shows the ensemble- 
averaged adjacent-gap ratio (r) as a function of the energy 
E/E^ for an isotropic speckle pattern of intensity Vq = E^, 
where Ea- is the correlation energy. The horizontal green line 
is the result for the Wigner-Dyson distribution {r)^jj, and 
the dashed black line the one for the Poisson distribution 
(r)p. The inset gives comparison between different system 
sizes. The vertical orange line indicates the position of the 
mobility edge Ec (the hatched rectangle represents the error- 
bar). The gray bar represents the value of Ec predicted in 
Ref. [ 4 ^ using transfer-matrix theory. 


and determine the energy threshold Ec which separates 
them, we compute the ratio of consecutive level spacings: 
r = min {i5„, (5n_i} / max{(5ji, The average over 

disorder realizations is known to be {r)yyjj — 0.5307 for 
the Wigner-Dyson distribution, and (r)p ~ 0.38629 for 
the Poisson distribution [d^. This statistical parameter 
was first introduced in Ref. in the context of 
many-body localization. In Fig. [3] we show the data 
corresponding to the disorder strength Vq = Ea-, where 
Ecr = Imcj^ is the correlation energy. We find that the 
ensemble average (r) changes rapidly from (r)p to 
as the energy increases. While in an infinite system 
one would have a sudden transition between the two 
statistics (with a third distribution exactly at Ec 0) , in 
a finite system we have a rapid but continuous crossover. 
For energies E < Eq, the data drift towards (r)p as 
the system size L increases since the localized wave 
functions are independent only for L —>■ 00 , while they 
drift in the opposite direction for E > Ec- The crossing 
of the curves corresponding to different system sizes 
indicates the critical energy (see inset of figure [3]). To 
pinpoint Ec we fit the data close to the transition with 
the scaling Ansatz (r) = g \{E — Ec) [5l[, where 

V is the critical exponent of the correlation length and 
g[x\ is the scaling function (universal up to a rescaling 
of the argument) which we Taylor expand up to second 
order. For the case of Fig. [31 from the best-fit analysis 
we obtain Ec = 0.576{10)E^, in quantitative agreement 
with the result of transfer-matrix theory from Ref. [i^ : 
Ec = 0.57Q{7)Ecr. For the critical exponent we obtain 
iz = 1.6(2), which is consistent with the prediction for 
the Anderson model: jz = 1.571(8) [5^. It is worth 
mentioning that in the energy regime E ^ Vq classical 
particles would be completely delocalized since the 
energy threshold tp for classical percolation in three- 
dimensional speckle patterns is extremely small, namely, 
tp ~ 10“^Vb [^. We consider also a red-detuned speckle 
field. Its distribution of intensities PrdiV) is the opposite 
of what corresponds to blue-detuned speckles, that is, 
Prdiy) = Phd{—V). The corresponding average value 
is (R(r)) = — Vq. At the disorder strength Vq = E^^, 
we obtain the mobility edge Ec = —0.81(4)iJcr, which 
(marginally) agrees with the result of transfer-matrix 
theory: Ec = —0.863(6) [i^. It is worth noticing 
that for blue-detuned speckles the mobility edge is well 
below the average intensity of the potential, while for 
red-detuned speckles it is instead above it. This strong 
asymmetry was already found in Ref. using transfer- 
matrix theory, but it was not captured by previous 
approximate calculations based on the self-consistent 
theory of localization. This means that predicting the 
position of the mobility edge requires quantitatively 
accurate methods. 

We now turn the discussion to concrete experimen¬ 
tal configurations. We first consider the setup where a 
single laser beam with wavelength A, propagating along 
the positive 2 axis (see Fig. [T|), is transmitted through 
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a diffusive plate and then focused onto the atomic cloud 
using a lens with focal length /. We assume the lens 
to be uniformly lit over a circular aperture of diameter 
D. The simplified numerical procedure employed before 
for isotropic models [U not apply in this case. 

The complex amplitude of the speckle field at the posi¬ 
tion r = {x, y, z) measured from the focal point can be 
computed using the Fresnel diffraction integral (d^ : 


^1 (r) = exp [i27r [z + /) /A] x 


f 

. {x-af + {y-pf 


\(.z + f) 


dad/3, (1) 


where ai(a,/3) is the complex field amplitude at the 
point 1 = (a,/3) just behind the focusing lens. The 
potential intensity is Vi(r) = \Ai (r)|^. Equation [3] was 
derived assuming paraxial approximation. Consistently, 
we will consider only small (positive) displacements from 
the focal point: x,y, z -C /, D. A convenient procedure 
to evaluate the Fresnel integral is to simulate the effect 
of a large number N of scattering centers randomly 
placed on the aperture . On the lens plane, one 

has: ai(a,/3) = exp(i(^„)(5(^) (1 - 1 „), where 

1„ = (a„,/3„) is the position of the nth scatterer, o„ 
is the modulus of the corresponding scattered wave, 
and (j)n is its phase, which has to be sampled from 
a uniform random distribution in the interval from 
—TT to TT. To simulate the effect of a uniform illu¬ 
mination, which is the case considered in this Rapid 
Communication, the moduli o„ have to be identically 
and independently distributed random variables, while 
the random positions of the scattering centers must 
fill the aperture circle uniformly [d^. Substituting 
the expression for oi (a, /3) in eq. O one obtains the 
complex field Ai (r) as the sum of wavelets propagating 
from the scattering center to the observation point. 
The field Ai (r) then has to be normalized to have 
the desired average intensity Vq- We verified that for 
N ^ 100 the resulting potential has the statistical 
properties of fully developed speckle patterns [i^. The 
intensities have the exponential probability distribution 
Rbd(R) defined above. The spatial correlation function 
is anisotropic, with cylindrical symmetry around the 
propagation axis z. If one takes two points aligned 
in the radial direction, the correlation function reads 
Frad(r) = [2Ji (r/cr) / (r/cr)]^ [i^, where the correlation 
length is fixed by the parameters of the optical appara¬ 
tus: tr = A// {Dtt). Instead, in the axial direction the 
correlation function is [d^: r^(r) = [siTi{r/az)/{r/cr, 
where the axial correlation length is Cz = 8X{f/D)'^/ tt. 
In current experimental implementations, the optical 
parameters are typically such that Cz > cr, meaning 
that the speckle grains are elongated along the beam 
propagation axis. For example, in the experiment of 
Ref. [l3| the anisotropy par ameter was ctz/o" ~ 6, while 
in the later experiment (55l| it was varied in the range 
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FIG. 4: (color online) Mobility edge Ec/Ea- as a function 
of the anisotropy parameter (Tz/ct. Uz and a are the axial 
and radial correlation lengths, respectively. The horizontal 
purple line indicates the disorder intensity Vo = Ea- The 
solid black curve is a guide to the eye (the dashed parts are 
an extrapolation). 


1 ^ ^T./cr < 10 by adjusting the aperture of the focusing 
length 

In order to investigate the effects due to the correlation 
anisotropy, we compute for varying values of the 
axial correlation length (Tz, considering both squeezed 
speckle grains {gzI<^ < 1) ^^nd elongated speckle grains 
[ozlo > I). In our computations the box shape is 
adapted to the disorder anisotropy, see [d^. The 
disorder intensity is set at Vq = = /i^/mcr^, defined 

using the (fixed) radial correlation length cr. We find 
that Ec monotonously decreases as we increase the 
anisotropy parameter Ozjo. In the quasi-isotropic case 
CTz/cr = I, the result agrees with the idealized isotropic 
model considered above, while it is approximately 50% 
larger for (Tz/ct = 1/9, and 15% lower for Uzjo = 6 
(see Fig. [3]). It is worth noticing that this dependence 
of Ec on the anisotropy parameter is not trivially 
related to the scaling of the average correlation energy 
E^ = defined from the geometric mean of 

the correlation lengths in the three spatial directions: 
tf = (crcrcTz)^^^. This suggests that the geometric mean 
tf is not the unique relevant length-scale, and that the 
structure of the spatial correlations plays a central role. 
We emphasize that in this Rapid Communication we 
are considering the speckle pattern created by a uniform 
aperture function. With different kinds of illumination 
(e.g, the Gaussian illumination [il[il), Ec might be 
somewhat different. 

While the reduction of Ec due to a large axial correlation 
length could be observed using currently available exper¬ 
imental setups, the increase of Ec is not easily accessible 
since the optical apparatuses do not permit us to create 
squeezed speckle grains. However, we can show that a 
similar increase in Ec is induced when two orthogonal 
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speckle patterns are superimposed. Explicitly, we 
numerically construct the potential due to the sum of 
two speckle patterns generated by laser beams with the 
same wavelength A. The first pattern propagates along 
z and the second along y (see Fig. [T|), and they interfere 
coherently, as is the case when the two laser beams 
have the same linear polarization. The total complex 
amplitude is then [13, Atot (r) = Ai (r) + A 2 (r). 
The complex amplitude of the second speckle pattern 
A 2 (r) can be computed using eq. [3] as described above, 
just exchanging the roles of the coordinates y and z 
in the right-hand side. For simplicity, we consider two 
patterns created with equal circular apertures, lit with 
the same (uniform) intensity, and focused using identical 
lenses. Thus the corresponding potentials |Ai(r)|^ and 
|A 2 (r)|^ have the same radial correlations lengths, which 
we set at tr = 0.75A. Their axial correlation lengths are 
extremely large, so that their variations along the respec¬ 
tive propagation axes are irrelevant. This configuration 
is inspired by the experimental setup of Ref. [42| . The 
potential E(r) = |4l(r)| corresponding to the coherent 
sum of two blue-detuned fields has the same exponential 
intensity distribution Pbd(R) as a single (blue-detuned) 
speckle pattern [d^. The structure of the spatial 
correlations of this total potential is instead much more 
intricate [i^. To describe it, it is convenient to consider 
the principal axes y' and z', obtained with a 45° rotation 
of the y and z axes around the x axis (see Fig. [T]). 
The correlation between two points aligned in parallel 
with the first principal axis z' is T^ (r) = r''ad(r/-\/2), 
meaning that the correlation length is dp = ^/2a (4^ . 
Moving in parallel with the second principal axis y', the 
potential is seen to oscillate rapidly due to the interfer¬ 
ence fringes and the corresponding correlation function 
is: (r) = [2Ji (r/cTp) / (r/cTp) cos (v^Trr/A)] The 

correlation function along the transverse axis x is instead 
the same as for a single speckle pattern: r“(r) = r''ad(r'). 
For our choice of parameters, the correlation function 
along the second principal axis T^ (r) touches zero 
four times before the first zero of the corresponding 
function along the transverse axis r^(7’), indicating 
the strong anisotropy of the disorder correlations. 
We consider again a potential with average intensity 
Vq = Efj. The corresponding mobility edge is found to 
be Ec = 0.67(l)£'cr, significantly higher than for the 
idealized isotropic disorder. This result is comparable 
with the one for a single speckle-pattern with squeezed 
axial correlation length (Jzja ~ 1/3 . We argue that 
this increase of E,. is induced by the rapid variations 
of the potential due to the interference fringes, which 
effectively reduce the spatial correlation length along 
the second principal axis. Experimentalists can easily 
modify the width of the interference fringes, either by 
changing the angle between the two beams or by using 
lasers with different wavelengths. Observing the increase 
of Ec is thus within experimental reach. 


We now turn the discussion to the comparison with 
the available experimental data. E^ was first measured 
in Ref. [l^ in the single-pattern configuration. The 
speckle grains were elongated, corresponding approxi¬ 
mately to the anisotropy parameter CT^/tT ~ 6 [^. In 
the regime of disorder strengths Vq ~ Fo-, the results 
were in the range 1.5Vb % E^ % 2Vo- These findings 
do not agree with our results for strongly elongated 
speckle grains: E^ ~ 0.5Vo- Most likely, the reason 
of this discrepancy traces back to the procedure used 
to extract the values of E^ from the measurement of 
the fraction of atoms that remain Anderson localized. 
In this derivation, the spectral function was approx¬ 
imated using the disorder-free value [s^ 1^ . This 
approximation is not reliable at the disorder strength 
necessary to observe Anderson localization. Notice also 
that in the experiment of Ref. [isj a Gaussian pupil 
function was employed. More recently, the mobility 
edge was measured in the configuration with two crossed 
speckle patterns created with approximately uniform 
apertures [l^ . For disorder intensities comparable to 
the correlation energy Vb ~ the mobility edge was 
found in the regime Ac ~ Vq- This result is significantly 
larger than the predictions for idealized isotropic models 
of the disorder and, in this sense, is consistent with our 
findings. However, it also overestimates our prediction 
for Vb = Act. This discrepancy is probably due to the 
fact that in the experiment the two interfering speckle 
patterns are not equivalent because they were created 
using slightly different apertures and lenses with differ¬ 
ent focal lengths. Also, the width of the experimental 
interference fringes is slightly smaller compared that in 
to our model. Furthermore, an exact modeling of the 
experiment of Ref. [l^ would require us to go beyond 
the paraxial approximation. All of these details of the 
experimental setup, once fully characterized, could be 
easily implemented in our formalism to compute Ac. 

In conclusion, we have studied the Anderson local¬ 
ization of matter waves exposed to optical speckles 
in the framework of quantum-chaos theory. We have 
shown that the structure of the spatial correlation of 
the disorder determines the position of the mobility 
edge, and we have described the effects induced by 
the correlation anisotropy in concrete experimental 
configurations, thus paving the way to a quantitative 
comparison between theory and experiment. This study 
sets the basis for future investigations of the effects due 
to interactions on the transport and on the coherence 
properties of disordered atomic gases and on the role 
played by the fractality of the critical wave functions 
close the mobility edge [^ . 
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Supplemental Material for 
Anderson localization of matter waves in 
quantum-chaos theory 


In the following supplemental material, we provide 
additional technical details about the numerical proce¬ 
dure we employ to determine the energy spectrum and 
the level-spacing statistics of isotropic and anisotropic 
speckle patterns. 

The real-space Hamiltonian of a quantum particle 
moving in a speckle pattern is given by: H = — -|- 

y(r), where h is the reduced Planck’s constant, m the 
particle’s mass, and V (r) is the external potential at the 
position r corresponding to the intensity of the optical 
speckle field. We consider a box with periodic boundary 
conditions and linear dimensions L^, Ly and L^., in the 
three directions l = x,y,z. 

It is convenient to represent the Hamiltonian operator 
in momentum space as a large hnite matrix: Hk.k' = 
7k,k' + Vic.k'; 2k,k' represents the kinetic energy op¬ 
erator, and Vk,k' the potential energy operator. The 
wavevectors form a discrete three-dimensional grid: k = 
{kx^ ky^ kz)^ with the three components k^ = where 
jt = — W/2,..., Ai/2 — 1, and the (even) number W 
determines the size of the grid in the 6 direction and, 
hence, the corresponding maximum wavevector. There¬ 
fore, when expanded in the square matrix format, the 
size of the Hamiltonian matrix iJk k' is N^ot x iVtot j where 
iVtot = NxNyNz. 

In this basis the kinetic energy operator is diagonal: 
2k,k' = where Sk,,k[ is the Kro- 

necker delta. The element 14,k' of the potential energy 
matrix can be computed as: 14,k' = Wk'-k, where Fk 
is the discrete Fourier transform of the speckle pattern 
F(r): 


Vka:kyk^ — ^ ^ ^ 

rx Ty 

exp [-i [k^Tx + kyTy -|- k^Tz)] I (2) 

here, =V{v = {rx,ry,rz)) is the value of the ex¬ 

ternal potential on the Ntot nodes of a regular lattice de¬ 
fined by Ty = LyifiJN^, with = 0,.., W — 1. Wavevec¬ 
tors differences are computed exploiting periodicity in 
wavevector space. 

The numerical algorithm we employ to generate the 
speckle patterns is described in the main text. It allows 
us to create both isotropic and anisotropic speckles. In 
the former case, the scattering centers (see main text) 
have to be placed on a fictitious spherical shell of diame¬ 
ter £), with uniform random distribution. Also, equation 
(1) of the main text has to be trivially modified, arriving 
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FIG. 5: Two-point spatial correlation functions of isotropic 
speckle patterns generated by scattering centers placed on a 
continuous spherical shell and on a discrete grid. The con¬ 
tinuous blue curve represents the analytical formula r‘'*°(r). 
The size of the cubic box is L = IGrrcr. 


at: 

di (r) = exp [i27r//A] x 
iXi 

On exp {i(j)n) exp 

n 

where 1„ = (a„,/3n,7„) is the position of the n-th scat¬ 
tering center on the (fictitious) spherical aperture. One 
obtains a speckle pattern with the isotropic correlation 
function r*®°(r) (see definition in the main text). It 
is worth mentioning that in order to generate isotropic 
speckle patterns one could instead employ the conven¬ 
tional numerical recipe described in Refs. [Sl|, [S^. 

For isotropic speckles, we employ cubic boxes with L = 
Lx = Ly = Lz (with periodic boundary conditions). The 
numerical recipe described here (and in the main text) 
creates a potential which does not necessarily satisfy pe¬ 
riodic boundary conditions. In principle, this could in¬ 
troduce distortions in the Fourier transform. However, a 
speckle pattern compatible with periodic boundary con¬ 
ditions is obtained if the position of the scattering centers 
is discretized according to: —?► nint (q;„/( 5) 5, where 

nint (•) is the function that returns the whole integer 
closest to the argument, <5 = A//L is the discretization 
step, and the analogous operation is applied to /?„ and 
7 „. In figure [SJ we compare the correlation functions 
numerically measured in speckle patterns generated with 
continuous and with discrete scatterers, against the ana¬ 
lytical formula r‘®°(r). The perfect agreement indicates 
that the discretization procedure does not alter the sta¬ 
tistical properties of the speckles. 

In figure [HI we compare the analyses of the energy-levels 
statistics of two speckle patterns (with average intensity 
Vq = E„) generated with continuous and with discrete 


FIG. 6: Analyses of the level-spacing statistics for isotropic 
speckle patterns generated by scattering centers placed on a 
continuous spherical shell and on a discrete grid, (r) is the 
ensemble-averaged adjacent-gap ratio, which is plotted as a 
function of the energy E/E^, where Ea is the correlation 
energy. The horizontal green line is the result for the Wigner- 
Dyson distribution (r)^^, the dashed black line the one for 
the Poisson distribution (r)p. 


scatterers. The excellent agreement cross-validates the 
two procedures, meaning that they are both suitable for 
our purposes. This is to be expected, since the linear size 
of the cubic box is L = 207r(T (this is a typical size we 
use), much larger than the correlation length of the disor¬ 
der, so that border effects play a minor role. Notice that 
the parameter cr characterizes the lengths scale of the dis¬ 
order spatial correlation, and hence, the typical speckle 
size. More quantitatively, the full width at half maxi¬ 
mum £c, defined by the condition T'^°{£c/2) = r‘®°(0)/2, 
is £c = 0.89707. The discretization procedure can be eas¬ 
ily adapted to have periodicity in anisotropic speckle pat¬ 
terns, provided the focal length / is much larger than the 
box size. 

A crucial step to guarantee the accuracy of our result is 
to test that the number of wavevectors (and, correspond¬ 
ingly, the maximum wavevector) is sufficient to have an 
accurate representation of the speckle pattern and of 
the orbitals. In figure [71 we compare the level-spacing 
statistics for an isotropic speckle pattern obtained with 
different numbers of wavevectors. For isotropic speckles, 
it is convenient to set Nx = Ny = = N. It is evident 

that we obtain convergence already for moderately 
large N ss 18. Consequently, we can afford to perform 
ensemble averages over a large number of disorder 
realizations in a suitable computational time. For 
example, the data corresponding to = 18 in figure [7] 
have been obtained by averaging approximately 3 x 10^ 
disorder realizations, requiring 72 hours on a 20-cores 
CPU, usin g th e PLASMA library for linear algebra com¬ 
putations (43| . The maximum wavevector number we 


|r - L 


ITT- 


A/ 


, (3) 

















9 



FIG. 7: Adjacent-gap ratio (r) as a function of the energy 
E/Ea, for an isotropic speckle pattern with intensity Vb = 
Ea- The different datasets correspond to different number of 
wavevectors N, thus to different levels of accuracy. 



FIG. 8: Adjacent-gap ratio (r) as a function of the energy 
E/Ea-, for an anisotropic speckle pattern in the single-beam 
configuration. The anisotropy parameter is = 2/9, 

while the disorder intensity is Vb = E^. The box sizes 
are Lx = Ly = T.Sttct and Lz = LxjS, meaning that the 
box shape has been only partially adapted to the disorder 
anisotropy. The different datasets correspond to different 
number of wavevectors Nz in the axial direction z, while 
in the radial directions the wave-vector number is fixed at 
Nx = Ny^ 12 . 


consider is iV = 40, allowing us to solve approximately 
13 disorder realizations in 24 hours. Notice that, when 
we increase the system size L, we proportionally increase 
N, so that the maximum wavevector in the Hamiltonian 
matrix remains fixed and, therefore, we maintain the 
same level of accuracy. 

In the case of anisotropic speckle patterns, it is conve¬ 
nient to (partially) adapt the shape of the box to the 



FIG. 9: Adjacent-gap ratio (r) as a function of the energy 
E/Ea, for the same anisotropic speckle pattern as in figure |8] 
Here, the different datasets correspond to different number of 
wavevectors Nx = Ny in the radial directions x and y, while 
in the axial direction z the wave-vector number is fixed at 
Nz = 14. 


disorder anisotropy, so that the system size can be set 
to be much larger than the disorder correlation length 
in each direction, without exceedingly increasing the 
matrix size. Also, the numbers of wavevectors in the 
three spatial directions N^, Ny, and have to be 
adjusted according to the corresponding linear system 
sizes. Lx, Ly, and L^, respectively, and also according 
to the correlation length in the corresponding direction. 
The shorter the correlation length, the larger Ni_ is 
required. In the single laser-beam configuration, the 
disorder correlations have axial symmetry around the 
beam propagation axis z; therefore we set Lx = Ly and 
Nx = Ny. As an illustrative example, we consider here 
an anisotropic pattern with strongly squeezed grains 
corresponding to ctz/ct = 2/9. The parameter Cz charac¬ 
terizes the correlation length in the axial direction; the 
full width at half maximum of the axial corresponding 
correlation function T^{z) (see definition in the main 
text) is ic = 0.897rcrz. For the radial correlation function 
one has ^c = 1.0297rcr. Notice that even in 
the case <j = az, the disorder is not perfectly isotropic. 
In our anisotropic example, we employ a box with an 
anisotropy which is similar to that of the disorder: 
Lz = Lx/?!. It is important to stress that the analysis 
of the effect due to the finite wavevectors number has 
to be performed both for the axial direction z and the 
radial directions x and y, separately. The former effect 
is analyzed in figure HI the latter in figure O Again, we 
observe convergence with moderately large grid sizes. 
We emphasize that the analysis of the level-spacing 
statistic can be used to determine the mobility edge 
both in the case of isotropic and anisotropic models; see, 
e.g.. Ref. [13 . In particular, it was found in Ref. [S5l| 
that the position of the mobility edge does not depend 
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FIG. 10: Main panel: ensemble-averaged adjacent-gap ra¬ 
tio (r) as a function of the energy E/Ea- for an anisotropic 
speckle pattern of intensity Vo — E^ and anisotropy parame¬ 
ter az/cr — 2/9. The different datasets correspond to differ¬ 
ent box sizes ~ Ly, with LzjLx = 1/3. Inset: compari¬ 
son between different system sizes. The vertical orange line 
indicates the position of the mobility edge E^ (the rectangle 
with pattern represents the error-bar). The continuous curves 
represent the scaling Ansatz g[x] (defined in the main text) 
expanded to second order. 


on the shape of the box, and that the level-spacing 
statistics converges in the thermodynamic limit to the 
Wigner-Dyson and to the Poisson distributions, in the 
delocalized and in the localized regimes, respectively. 
Instead, the level-spacing statistics exactly at the critical 
point E = Ec (which is system-size independent) was 
found to depend on the box shape. This does not 
invalidate our finite-size scaling procedure, but likely 
implies that the amplitude of the parameter (r) at 
Ec for different speckle anisotropies is not universal. 
Results analogous to those of Ref. [S^ were obtained 
in Refs. [H, [S^, where the effect due to the choice of 
boundary conditions (e.g., periodic boundary conditions 
vs. hard-wall or Dirichlet boundary conditions) was an¬ 
alyzed. Again, it was found that the position of mobility 
edge is not affected by the choice of boundary conditions, 
while the level-spacing distribution at the critical point 
is. Figure [10] shows the analysis of the level-spacing 
statistics for our example of anisotropic speckle, and 


the finite-size scaling analysis (see inset) employed to 
locate the mobility edge. For completeness, we have also 
analyzed the potential effect on Ec due to the box shape 
by repeating calculations with different LjLz ratios. 
Consistently with the findings of Refs. [H, [H, [S^, we 
also find that there is no systematic bias by changing 
the box shape; however, larger errorbars are obtained 
if the shape is not the optimal one, since one has to 
employ larger matrix sizes. 

In the configuration with two superimposed speckle 
patterns, the disorder has the intricate anisotropic 
spatial correlation structure described in the main text; 
see figure (1) [panel (c)] and figure (2) in the main text. 
In this case, we use a cubic box with = Ly = L^, 
but we rotate the speckle pattern such that second 
principal axis y' (see panel (d) of figure (1) in the main 
text) is aligned with one of the sides of the box, e.g., 
along the 2 direction. Then, we increase the wavevector 
number corresponding to this spatial direction, since the 
disorder has rapid oscillations and a shorter effective 
correlation length due to the interference fringes. For 
the two-pattern configuration described in the main 
text, we find it is necessary to use a wavevector number 
Nz « 3N,. 


Finally, it is worth comparing the efficiency of our pro¬ 
cedure to determine E^ which is based on the analysis 
of the level-spacing statistics within quantum-chaos the¬ 
ory, with the one of the transfer-matrix theory employed 
in Ref. [S^ in the case of isotropic speckle patterns. In 
the main text we show the quantitative comparison be¬ 
tween the two theories both for blue-detuned and red- 
detuned isotropic speckle fields. The results of transfer- 
matrix theory have somewhat smaller errorbars, in par¬ 
ticular for red-detuned speckles, since in our formalism 
they require a larger wavevector number, and hence, al¬ 
low us to consider fewer disorder realizations. However, 
we have shown that quantum-chaos theory provides us 
with effective and ffexible tools to address disorder pat¬ 
terns with intricate correlation structures and to adapt 
the shape of the sample to the disorder anisotropy. More 
importantly, our exact diagonalization study allows us 
to disclose interesting properties of the energy spectrum 
of optical speckles, thus creating a strong link between 
random matrix theory and ultracold atomic gases. 
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